{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a571e3ac",
   "metadata": {},
   "source": [
    "# Build cutouts for modeling\n",
    "\n",
    "\n",
    "This notebook contains 5 cells:\n",
    "1. Imports + defining main directory where data is stored\n",
    "2. utility functions that handle coordinte transformation and file handling\n",
    "3. additoinal function for running deepCR on case by case basis (might be useful for some ACS imaging)\n",
    "4. Main function that builds cutout files and writs as HDF5 and FITS formats\n",
    "5. Read proposal specific csv file and build cutout list\n",
    "6. loop over cutout list to build cutout files for each target/filter\n",
    "\n",
    "\n",
    "\n",
    "Last ran 2025.10.26 by C. Watson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "79e4e5e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% Imports\n",
    "import os, re, csv, json, glob, shutil\n",
    "from pathlib import Path\n",
    "from typing import Optional, Tuple\n",
    "import numpy as np\n",
    "import h5py\n",
    "from astropy.io import fits\n",
    "from astropy.wcs import WCS\n",
    "from astropy.nddata import Cutout2D\n",
    "from astropy.stats import sigma_clipped_stats\n",
    "from astropy.coordinates import SkyCoord\n",
    "import pandas as pd\n",
    "from astropy import units as u\n",
    "from deepCR import deepCR # use if properly installed\n",
    "from astropy.visualization import simple_norm\n",
    "\n",
    "Main_dir = '/Volumes/AGEL/HST_data'  # Directory where data is stored\n",
    "\n",
    "def final_irimg_name(target, propno, band, camera):\n",
    "    final_ir_image = f\"{Main_dir}/{target}/HST/{propno}_{band}/{target}_{band}_{camera}_irdrz_L1.fits\"\n",
    "    return final_ir_image\n",
    "\n",
    "def final_uvimg_name(target, propno, band, camera):\n",
    "    final_uv_image = f\"{Main_dir}/{target}/HST/{propno}_{band}/{target}_{band}_{camera}_uvdrz_L1.fits\"\n",
    "    return final_uv_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3687376",
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Utilities\n",
    "def _load_hdf5_2d(path, dataset=None):\n",
    "    \"\"\"Return a 2D numpy array from an HDF5 file.\n",
    "    If dataset is None, try common keys then fall back to the largest 2D dataset.\n",
    "    \"\"\"\n",
    "    with h5py.File(path, \"r\") as f:\n",
    "        if dataset:\n",
    "            if dataset not in f:\n",
    "                raise KeyError(f\"Dataset '{dataset}' not found in {path}\")\n",
    "            arr = f[dataset][...]\n",
    "            if arr.ndim != 2:\n",
    "                raise ValueError(f\"Dataset '{dataset}' is {arr.ndim}D, expected 2D\")\n",
    "            return np.asarray(arr, dtype=\"float64\")\n",
    "\n",
    "        # try common names first\n",
    "        for k in (\"kernel\", \"psf\", \"psf_kernel\", \"data\", \"image\"):\n",
    "            if k in f and isinstance(f[k], h5py.Dataset) and f[k].ndim == 2:\n",
    "                return np.asarray(f[k][...], dtype=\"float64\")\n",
    "\n",
    "        # fallback: find largest 2D dataset\n",
    "        best = None; best_n = -1\n",
    "        def walk(g):\n",
    "            nonlocal best, best_n\n",
    "            for name, obj in g.items():\n",
    "                if isinstance(obj, h5py.Dataset) and obj.ndim == 2:\n",
    "                    n = int(np.prod(obj.shape))\n",
    "                    if n > best_n:\n",
    "                        best, best_n = obj, n\n",
    "                elif isinstance(obj, h5py.Group):\n",
    "                    walk(obj)\n",
    "        walk(f)\n",
    "        if best is None:\n",
    "            raise ValueError(f\"No 2D dataset found in {path}\")\n",
    "        return np.asarray(best[...], dtype=\"float64\")\n",
    "\n",
    "def _load_fits_2d(path):\n",
    "    with fits.open(path, memmap=False) as hdul:\n",
    "        for h in hdul:\n",
    "            if h.data is None:\n",
    "                continue\n",
    "            arr = np.asarray(h.data)\n",
    "            if arr.ndim == 2:\n",
    "                return arr.astype(\"float64\")\n",
    "    raise ValueError(f\"{path} had no 2D image HDU\")\n",
    "\n",
    "def _load_psf_any(psf_path, dataset=None, normalize=False, ensure_odd=False, center_peak=False):\n",
    "    \"\"\"Load PSF from FITS or HDF5; optional normalize/shape/centering.\"\"\"\n",
    "    if psf_path is None:\n",
    "        return None\n",
    "    p = str(psf_path).lower()\n",
    "    if p.endswith((\".fits\", \".fit\", \".fz\")):\n",
    "        psf = _load_fits_2d(psf_path)\n",
    "    elif p.endswith((\".h5\", \".hdf5\")):\n",
    "        psf = _load_hdf5_2d(psf_path, dataset=dataset)\n",
    "    else:\n",
    "        raise ValueError(f\"Unsupported PSF format: {psf_path}\")\n",
    "\n",
    "    # optional fixes\n",
    "    if ensure_odd:\n",
    "        ny, nx = psf.shape\n",
    "        if ny % 2 == 0 or nx % 2 == 0:\n",
    "            psf = np.pad(psf, ((0, ny % 2 == 0), (0, nx % 2 == 0)), mode=\"constant\")\n",
    "\n",
    "    if center_peak:\n",
    "        try:\n",
    "            from scipy.ndimage import fourier_shift\n",
    "            yy, xx = np.unravel_index(np.nanargmax(psf), psf.shape)\n",
    "            cy, cx = (np.array(psf.shape) - 1) / 2.0\n",
    "            shift = (cy - yy, cx - xx)\n",
    "            psf = np.fft.ifftn(fourier_shift(np.fft.fftn(psf), shift)).real\n",
    "        except Exception:\n",
    "            pass\n",
    "\n",
    "    if normalize:\n",
    "        s = np.nansum(psf)\n",
    "        if not np.isfinite(s) or s == 0:\n",
    "            raise ValueError(\"PSF sum is zero/NaN; cannot normalize.\")\n",
    "        psf = psf / s\n",
    "\n",
    "    return psf\n",
    "\n",
    "def _load_weight_any(weight_path, dataset=None):\n",
    "    \"\"\"Load a 2D weight / ivar / rms map from FITS or HDF5.\"\"\"\n",
    "    if weight_path is None:\n",
    "        return None\n",
    "    p = str(weight_path).lower()\n",
    "    if p.endswith((\".fits\", \".fit\", \".fz\")):\n",
    "        return _load_fits_2d(weight_path)\n",
    "    elif p.endswith((\".h5\", \".hdf5\")):\n",
    "        return _load_hdf5_2d(weight_path, dataset=dataset)\n",
    "    else:\n",
    "        raise ValueError(f\"Unsupported weight format: {weight_path}\")\n",
    "\n",
    "def _read_cd_or_pc(header, force_pc: bool=False) -> Tuple[float,float,float,float]:\n",
    "    \"\"\"Return a linear WCS 2x2 in deg/pix from CD or PC*CDELT.\"\"\"\n",
    "    def has_cd(h): return all(k in h for k in (\"CD1_1\",\"CD1_2\",\"CD2_1\",\"CD2_2\"))\n",
    "    def has_pc(h): return all(k in h for k in (\"PC1_1\",\"PC1_2\",\"PC2_1\",\"PC2_2\")) and (\"CDELT1\" in h and \"CDELT2\" in h)\n",
    "    if (not force_pc) and has_cd(header):\n",
    "        return float(header[\"CD1_1\"]), float(header[\"CD1_2\"]), float(header[\"CD2_1\"]), float(header[\"CD2_2\"])\n",
    "    elif has_pc(header):\n",
    "        CDELT1, CDELT2 = float(header[\"CDELT1\"]), float(header[\"CDELT2\"])\n",
    "        return (CDELT1*float(header[\"PC1_1\"]), CDELT1*float(header[\"PC1_2\"]),\n",
    "                CDELT2*float(header[\"PC2_1\"]), CDELT2*float(header[\"PC2_2\"]))\n",
    "    else:\n",
    "        raise ValueError(\"No CD or PC*CDELT found in header.\")\n",
    "\n",
    "def _build_transform_pix2angle_arcsec(CD1_1, CD1_2, CD2_1, CD2_2) -> np.ndarray:\n",
    "    return np.array([[CD1_1, CD1_2],[CD2_1, CD2_2]], float) * 3600.0  # arcsec/pix\n",
    "\n",
    "def _centered_frame_ra_dec_at_xy0(nx: int, ny: int, A_arcsec: np.ndarray) -> Tuple[float,float]:\n",
    "    \"\"\"Origin at cutout center → center maps to (0,0) arcsec.\"\"\"\n",
    "    x_c, y_c = int(nx/2), int(ny/2)\n",
    "    dra, ddec = A_arcsec.dot(np.array([x_c, y_c], float))\n",
    "    return -float(dra), -float(ddec)\n",
    "\n",
    "def _estimate_background_rms(data: np.ndarray, sigma: float=3.0, maxiters: int=5) -> float:\n",
    "    _, med, std = sigma_clipped_stats(data, sigma=sigma, maxiters=maxiters, grow=False)\n",
    "    return float(std)\n",
    "\n",
    "def _load_first_2d(path: Optional[str]) -> Optional[np.ndarray]:\n",
    "    if not path: return None\n",
    "    with fits.open(path, memmap=False) as hdul:\n",
    "        for h in hdul:\n",
    "            if h.data is None: continue\n",
    "            arr = np.asarray(h.data)\n",
    "            if arr.ndim == 2:\n",
    "                return arr.astype(\"float64\")\n",
    "    raise ValueError(f\"{path} had no 2D image HDU\")\n",
    "\n",
    "def _resolve_center_pixels(header, wcs: WCS, center_ra=None, center_dec=None, center_x=None, center_y=None):\n",
    "    if center_x is not None and center_y is not None: return float(center_x), float(center_y)\n",
    "    if center_ra is not None and center_dec is not None:\n",
    "        x, y = wcs.world_to_pixel_values(center_ra, center_dec)\n",
    "        return float(x), float(y)\n",
    "    nx, ny = int(header.get(\"NAXIS1\")), int(header.get(\"NAXIS2\"))\n",
    "    return float(nx/2.0), float(ny/2.0)\n",
    "\n",
    "def _size_pixels_from_arcsec(size_arcsec: float, A_arcsec: np.ndarray) -> int:\n",
    "    eff_scale = float(np.sqrt(np.linalg.det(A_arcsec.T @ A_arcsec)) / np.sqrt(2.0))\n",
    "    if not np.isfinite(eff_scale) or eff_scale <= 0:\n",
    "        eff_scale = float(np.median(np.abs(A_arcsec)))\n",
    "    n_pix = int(np.round(size_arcsec / eff_scale))\n",
    "    return max(8, n_pix)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cc103e29",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Core function: write lenstronomy-ready HDF5 (+ optional FITS+PNG)\n",
    "def build_cutout_h5(\n",
    "    fits_path: str, output_h5: str,\n",
    "    size_pixels: int|None=None, size_arcsec: float|None=None,\n",
    "    center_ra: float|None=None, center_dec: float|None=None,\n",
    "    center_x: float|None=None, center_y: float|None=None,\n",
    "    psf_path: str|None=None, weight_path: str|None=None,\n",
    "    psf_h5_dataset: str|None=None, weight_h5_dataset: str|None=None,\n",
    "    psf_normalize: bool=False, psf_ensure_odd: bool=False, psf_center_peak: bool=False,\n",
    "    background_rms: float|None=None, exposure_time: float|None=None,\n",
    "    force_pc: bool=False, overwrite: bool=False,\n",
    "    save_fits_with_header: bool=False, thumbnail_png: str|None=None,\n",
    "    z_src: float|None=None, z_def: float|None=None, \n",
    "    object_name: str |None=None, propno: float|None=None, band: str|None=None, camera: str|None=None\n",
    "):\n",
    "    if (size_pixels is None) and (size_arcsec is None):\n",
    "        raise ValueError(\"Specify either size_pixels or size_arcsec.\")\n",
    "    if (center_ra is not None) ^ (center_dec is not None):\n",
    "        raise ValueError(\"Provide both center_ra and center_dec, or neither.\")\n",
    "    if (center_x is not None) ^ (center_y is not None):\n",
    "        raise ValueError(\"Provide both center_x and center_y, or neither.\")\n",
    "\n",
    "    with fits.open(fits_path, memmap=False) as hdul:\n",
    "        hdu = next(h for h in hdul if (h.data is not None and h.data.ndim==2))\n",
    "        header = hdu.header\n",
    "        img = np.asarray(hdu.data).astype(\"float64\")\n",
    "        w = WCS(header)\n",
    "        CD1_1, CD1_2, CD2_1, CD2_2 = _read_cd_or_pc(header, force_pc=force_pc)\n",
    "        A_arcsec = _build_transform_pix2angle_arcsec(CD1_1, CD1_2, CD2_1, CD2_2)\n",
    "        cx, cy = _resolve_center_pixels(header, w, center_ra, center_dec, center_x, center_y)\n",
    "        if size_pixels is None:\n",
    "            size_pixels = _size_pixels_from_arcsec(size_arcsec, A_arcsec)\n",
    "        cut = Cutout2D(img, (cx, cy), size=(int(size_pixels), int(size_pixels)), wcs=w, mode=\"trim\", copy=True)\n",
    "        data = cut.data.astype(\"float64\")\n",
    "        ny, nx = data.shape\n",
    "        ra0, dec0 = _centered_frame_ra_dec_at_xy0(nx, ny, A_arcsec)\n",
    "        bkg = float(background_rms) if background_rms is not None else _estimate_background_rms(data)\n",
    "        exptime = float(exposure_time) if exposure_time is not None else float(header.get(\"EXPTIME\", 1.0))\n",
    "        # kernel_psf = _load_first_2d(psf_path)\n",
    "        kernel_psf = _load_psf_any(psf_path,\n",
    "                               dataset=psf_h5_dataset,\n",
    "                               normalize=psf_normalize,\n",
    "                               ensure_odd=psf_ensure_odd,\n",
    "                               center_peak=psf_center_peak)\n",
    "        # weight_map = _load_first_2d(weight_path)\n",
    "        weight_map = _load_weight_any(weight_path, dataset=weight_h5_dataset)\n",
    "\n",
    "\n",
    "    out_path = Path(output_h5)\n",
    "    if (not overwrite) and out_path.exists():\n",
    "        raise FileExistsError(f\"{output_h5} exists. Set overwrite=True.\")\n",
    "    out_path.parent.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "    with h5py.File(out_path, \"w\") as f:\n",
    "        f.create_dataset(\"image_data\", data=data, dtype=\"float64\")\n",
    "        f.create_dataset(\"background_rms\", data=np.array(bkg, dtype=\"float64\"))\n",
    "        f.create_dataset(\"exposure_time\", data=np.array(exptime, dtype=\"float64\"))\n",
    "        f.create_dataset(\"ra_at_xy_0\", data=np.array(ra0, dtype=\"float64\"))\n",
    "        f.create_dataset(\"dec_at_xy_0\", data=np.array(dec0, dtype=\"float64\"))\n",
    "        f.create_dataset(\"transform_pix2angle\", data=A_arcsec, dtype=\"float64\")\n",
    "        f.create_dataset(\"z_spec_SRC\", data=np.array(z_src if z_src is not None else -1.0, dtype=\"float64\"))\n",
    "        f.create_dataset(\"z_spec_DE\", data=np.array(z_def if z_def is not None else -1.0, dtype=\"float64\"))\n",
    "        if kernel_psf is not None:\n",
    "            f.create_dataset(\"kernel_psf\", data=kernel_psf, dtype=\"float64\")\n",
    "        if weight_map is not None:\n",
    "            f.create_dataset(\"weight_map\", data=weight_map, dtype=\"float64\")\n",
    "        meta = {\n",
    "            \"source_fits\": os.path.abspath(fits_path),\n",
    "            \"center_input\": {\"ra_deg\": center_ra, \"dec_deg\": center_dec, \"x_pix\": center_x, \"y_pix\": center_y},\n",
    "            \"size_pixels\": int(size_pixels),\n",
    "            \"CD_matrix_deg_per_pix\": [CD1_1, CD1_2, CD2_1, CD2_2],\n",
    "            \"transform_pix2angle_arcsec\": A_arcsec.tolist(),\n",
    "            \"background_rms_method\": \"sigma_clipped_stats\" if background_rms is None else \"user_provided\",\n",
    "            \"exptime_source\": \"FITS[EXPTIME]\" if exposure_time is None else \"user_provided\",\n",
    "            \"psf_attached\": bool(kernel_psf is not None),\n",
    "            \"weight_attached\": bool(weight_map is not None),\n",
    "            \"z_spec_SRC\": float(z_src) if z_src is not None else None,\n",
    "            \"z_spec_DE\": float(z_def) if z_def is not None else None\n",
    "        }\n",
    "        f.create_dataset(\"meta\", data=json.dumps(meta))\n",
    "\n",
    "    if save_fits_with_header:\n",
    "        # regenerate a clean WCS header for the cutout FITS and add LNS_* keys\n",
    "        with fits.open(fits_path, memmap=False) as hduL1:\n",
    "            h2 = next(h for h in hduL1 if (h.data is not None and h.data.ndim==2))\n",
    "            w2 = WCS(h2.header)\n",
    "            cut2 = Cutout2D(h2.data, (cx, cy), size=(int(size_pixels), int(size_pixels)), wcs=w2, mode=\"trim\", copy=True)\n",
    "        hdr = cut2.wcs.to_header()\n",
    "        hdr[\"LNS_RA0\"] = (ra0,  \"lenstronomy ra_at_xy_0 (arcsec)\")\n",
    "        hdr[\"LNS_DE0\"] = (dec0, \"lenstronomy dec_at_xy_0 (arcsec)\")\n",
    "        hdr[\"LNS_T11\"] = (A_arcsec[0,0], \"transform_pix2angle[0,0] (arcsec/pix)\")\n",
    "        hdr[\"LNS_T12\"] = (A_arcsec[0,1], \"transform_pix2angle[0,1] (arcsec/pix)\")\n",
    "        hdr[\"LNS_T21\"] = (A_arcsec[1,0], \"transform_pix2angle[1,0] (arcsec/pix)\")\n",
    "        hdr[\"LNS_T22\"] = (A_arcsec[1,1], \"transform_pix2angle[1,1] (arcsec/pix)\")\n",
    "        hdr[\"EXPTIME\"] = (exptime, \"exposure time (s)\")\n",
    "        if z_src is not None:\n",
    "            hdr[\"Z_SRC\"] = (float(z_src), 'source z_spec (AGEL Airtable)')\n",
    "        if z_def is not None:\n",
    "            hdr[\"Z_DEF\"] = (float(z_def), 'deflector z_spec (AGEL Airtable)')\n",
    "        fits_name = str(out_path.with_suffix(\".fits\"))\n",
    "        fits.PrimaryHDU(data=data, header=hdr).writeto(fits_name, overwrite=True)\n",
    "\n",
    "    return str(out_path)\n",
    "\n",
    "# Catalog runner\n",
    "def _float_or_none(v):\n",
    "    try:\n",
    "        if v is None: return None\n",
    "        s = str(v).strip()\n",
    "        if s == \"\" or s.lower() in {\"none\",\"nan\",\"null\"}: return None\n",
    "        return float(s)\n",
    "    except Exception:\n",
    "        return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bb134ff7",
   "metadata": {},
   "outputs": [],
   "source": [
    "## read in csv for specific proposal id and build cutout list\n",
    "proposal_id_no = '17307'\n",
    "cutout_lens = pd.read_csv(f'{Main_dir}/{proposal_id_no}_targets.csv')\n",
    "filters = ['F606W']#['F140W', 'F200LP']\n",
    "camera = 'ACS' \n",
    "\n",
    "cutout=[]\n",
    "for index, row in cutout_lens.iterrows():\n",
    "    \n",
    "    for band in filters:\n",
    "        target_name = row['objname']\n",
    "        if target_name[-1] == 'B':\n",
    "            continue  # skip B images for now\n",
    "        if target_name == 'AGEL085331+232155A' or target_name == 'AGEL014556+040229A' or target_name == 'AGEL015009-030438A':\n",
    "            continue  # skip this one for now - has bad coords/exposures\n",
    "        fits_file = glob.glob(f'{Main_dir}/{target_name}/HST/{proposal_id_no}_{band}/*_sci.fits')\n",
    "        if len(fits_file) == 0:\n",
    "            print(f'No fits file found for {target_name} in band {band}')\n",
    "            continue\n",
    "\n",
    "        if band == 'F140W':\n",
    "            pix_scale = 0.08\n",
    "        elif band == 'F200LP' or band == 'F606W':\n",
    "            pix_scale = 0.05\n",
    "        fits_file = fits_file[0]\n",
    "        hdul = fits.open(fits_file)\n",
    "\n",
    "        ra_deg = row['RAJ2000'] \n",
    "        dec_deg = row['DECJ2000'] \n",
    "        z_src = row['z_spec_SRC_Spectral_Observations_Tally'] if not np.isnan(row['z_spec_SRC_Spectral_Observations_Tally']) else row['z_source (from DR2 Redshifts)'] if not np.isnan(row['z_source (from DR2 Redshifts)']) else row['z_spec_SRC'] if not np.isnan(row['z_spec_SRC']) else None\n",
    "        z_def = row['z_spec_DE_Spectral_Observations_Tally'] if not np.isnan(row['z_spec_DE_Spectral_Observations_Tally']) else row['z_deflector (from DR2 Redshifts)'] if not np.isnan(row['z_deflector (from DR2 Redshifts)']) else row['z_spec_DE'] if not np.isnan(row['z_spec_DE']) else None\n",
    "        \n",
    "            \n",
    "        position=SkyCoord(ra_deg*u.deg,dec_deg*u.deg)\n",
    "        ra_tup = position.ra.hms\n",
    "        dec_tup = position.dec.dms\n",
    "        ra_list = [int(float(position.ra.hms[i])) if len(str(int(float(position.ra.hms[i]))))==2 else '0{}'.format(int(float(position.ra.hms[i]))) for i in range(len(ra_tup)-1)]\n",
    "        ra_list.append(float(position.ra.hms[2]))\n",
    "        dec_list = [int(float(str(position.dec.dms[i]).strip('-+'))) if len(str(int(float(str(position.dec.dms[i]).strip('-+')))))==2 else '0{}'.format(int(float(str(position.dec.dms[i]).strip('-+')))) for i in range(len(ra_tup)-1)]\n",
    "        dec_list.append(float(str(position.dec.dms[2]).strip('-+')))\n",
    "        ra = '{} {} {}'.format(ra_list[0], ra_list[1], ra_list[2])\n",
    "        if str(dec_tup[0])[0] == '-':\n",
    "            dec = '-{} {} {}'.format(dec_list[0], dec_list[1], dec_list[2])\n",
    "        else:\n",
    "            dec = '+{} {} {}'.format(dec_list[0], dec_list[1], dec_list[2])\n",
    "        cutout.append([row['objname'],\n",
    "                    band,\n",
    "                    camera,\n",
    "                    ra,\n",
    "                    dec,\n",
    "                    pix_scale,\n",
    "                    z_src,\n",
    "                    z_def])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "242a2dcd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done for AGEL051603-220847A F606W\n"
     ]
    }
   ],
   "source": [
    "# Loop over cutout list and build cutouts\n",
    "for i in cutout:\n",
    "    object_name = i[0]\n",
    "    band = i[1]\n",
    "    camera = i[2]\n",
    "    center_ra, center_dec = i[3], i[4]\n",
    "    fits_path = f\"{Main_dir}/{object_name}/HST/{proposal_id_no}_{band}/{object_name}_{band}_{camera}_{'drz' if band=='F140W' else 'drc'}_sci.fits\"\n",
    "    output_h5 = f\"{Main_dir}/{object_name}/HST/{proposal_id_no}_{band}/{object_name}_{band}_{camera}_cutout_L3.h5\"\n",
    "    sc = SkyCoord(center_ra, center_dec, unit=(u.hourangle, u.deg))\n",
    "    center_ra, center_dec = sc.ra.deg, sc.dec.deg\n",
    "    center_x, center_y = None, None\n",
    "\n",
    "    size_arcsec = 25.0  # desired cutout size in arcsec\n",
    "\n",
    "    size_pixels = size_arcsec / i[5]  # convert to pixels \n",
    "    psf_path = f\"/Volumes/AGEL/lens_processing/psf_model_{band}.h5\"\n",
    "    weight_path = f'{Main_dir}/{object_name}/HST/{proposal_id_no}_{band}/{object_name}_{band}_{camera}_{'drz' if band=='F140W' else 'drc'}_wht.fits'\n",
    "    build_cutout_h5(\n",
    "        fits_path, output_h5,\n",
    "        size_pixels=size_pixels, size_arcsec=size_arcsec,\n",
    "        center_ra=center_ra, center_dec=center_dec, center_x=center_x, center_y=center_y,\n",
    "        psf_path=psf_path, weight_path=weight_path,\n",
    "        overwrite=True, save_fits_with_header=True, thumbnail_png=None,\n",
    "        z_src=i[6], z_def=i[7], object_name = object_name, propno=proposal_id_no, band=band, camera=camera\n",
    "    )\n",
    "    print(f'Done for {object_name} {band}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fe446f7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "stenv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
